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ABSTRACT 


This  paper  reports  on  continued  research  toward  establishing  a  consistent  modeling  framework  for  calculating 
nuclear  explosions  in  earth  materials.  The  model  must  be  consistent  with  observed  phenomena  in  the  near-field  by 
correctly  1)  calculating  the  resulting  explosive  cavity  radius  for  a  given  yield  and  depth  of  burial,  2)  accounting  for 
the  correct  energy  deposition  by  partitioning  it  into  internal  (heat  and  plastic  strain)  and  kinetic  (e.g.  radiated 
seismic)  energy,  3)  predicting  the  free-field  displacement/velocities  waveforms  and  4)  predicting  the  measured 
attenuation  of  the  free-field  peak  velocity  with  distance.  The  model  developed  in  the  last  year  satisfies  all  of  these 
criteria  and  has  been  exercised  in  the  investigation  of  the  2009  North  Korean  nuclear  test.  The  main  findings 
reported  in  this  paper  are:  a)  the  extension  of  the  developed  model  to  the  analysis  of  scaled  depth  of  burial  and  free 
surface  effects  in  2D  Earth  structure,  and  b)  the  improvement  of  the  computational  equation  of  state  (EOS)  for 
granite/granodiorite  and  some  examples  of  the  models  self-consistency.  Our  study  focuses  on  the  North  Korean  test 
site  and  the  May  2009  test.  When  compared  to  the  Denny  and  Johnson  (1991)  and  to  the  Heard  and  Ackerman 
(1967)  cavity  radius  scaling  models,  the  results  presented  in  this  paper  show  a  clear  preference  to  the  statistical 
model  developed  by  Denny  and  Johnson.  In  addition,  comparative  work  between  Patton  (201 1)  and  the  model 
developed  under  this  project  provides  a  lower  limit  to  the  yield  and  depth  of  burial  for  the  2009  North  Korean  test.  A 
series  of  sensitivity  analysis  comprising  the  variation  of  key  material  properties  and  the  incorporation  of  topography 
is  being  produced  at  the  time  of  writing  this  paper  and  will  be  presented  at  the  conference.  This  extended  analysis 
will  provide  additional  bounds  on  the  uncertainty  for  these  estimates. 

A  strong  motion  hydrodynamics  code  was  used  to  investigate  the  dependence  of  cavity  dynamics  and  final  cavity 
radius  on  the  main  material  properties  (Young’s  modulus,  shear  modulus,  porosity,  etc.).  The  material  model 
developed  was  obtained  by  taking  the  Piledriver  and  the  Hardhat  nuclear  test  events  as  the  main  design  references. 
The  following  features  of  the  problem  were  identified  when  developing  the  material  model:  velocity  profiles  at 
given  stations  (near  field),  source  modeling  alternatives  (iron  pill,  ideal  gas,  Hydses/SESAME),  energy  partition 
after  the  shot,  peak  velocity  and  peak  displacement  attenuation  profiles  and  final  cavity  size  as  a  function  of  the 
depth  of  burial.  Previous  attempts  made  with  existing  material  models  failed  to  comply  with  one  or  more  of  these 
features.  The  Tillotson  type  of  equation  of  state  combined  with  a  shear  plastic  strain-dependent  strength  model  was 
implemented  and  used  to  observe  surface  effects  from  various  scaled  depths  of  burial  for  a  nuclear  explosion. 

The  developed  material  model  was  used  in  a  set  of  2D  axially  symmetric  simulations  with  a  flat  free  surface,  i.e.,  no 
topography.  Given  the  best  estimates  for  the  material  parameters  and  the  fact  that  there  was  no  evidence  of  surface 
spall  at  the  North  Korean  test  site  after  the  2009  nuclear  test,  the  calculations  place  a  minimum  yield  and  depth  of 
burial  of  5.7  kilotons  and  375  meters  for  a  uniform  source  region.  Further  refinement  of  these  numbers  will  only  be 
possible  by  introducing  a  more  realistic  topography  profile  corresponding  to  the  North  Korean  test  site. 
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OBJECTIVES 


Source  characterization  is  the  focus  of  this  project.  Understanding  the  physical  basis  of  seismic  wave  generation  is 
recognized  as  the  key  to  advancing  our  ability  to  characterize  the  explosion  source. 

The  first  main  objective  of  the  current  research  is  to  undertake  a  series  of  computer  simulations  to  study  the 
evolution  of  the  cavity  radius  under  different  conditions.  This  phase  of  computer  simulations  will  employ  LANL 
strong  motion  codes  (CASH)  for  uniform  source  region  structures  to  investigate  the  dependence  of  cavity  dynamics 
and  final  radius  on  material  properties,  such  as  Young’s  and  shear  moduli,  gas  porosity,  overburden  (pgh),  regional 
stresses,  etc.  The  calculations  were  done  using  a  material  having  some  similarity  to  the  granodiorite  that  exists  at  the 
Democratic  People’s  Republic  of  Korea  (DPRK)  test  site. 

The  second  main  objective  of  the  current  work  is  to  perform  2D  axially  symmetric  and  3D  computer  simulations 
that  will  introduce  sub-surface  heterogeneity  and  free  surface  topography  into  the  models.  Here  the  focus  will  be  to 
identify  the  impact  of  certain  idealized  source  region  structures  on  the  dynamics  of  cavity  growth  and  on  the 
radiated  seismic  energy.  These  structures  will  be  introduced  into  models  of  the  source  region  in  controlled  ways  in 
order  to  isolate  their  effects  on  the  close-in  phenomenology  and  to  identify  features  in  the  cavity  growth  and  seismic 
radiation  associated  with  those  effects. 

RESEARCH  ACCOMPLISHED 


Introduction.  A  number  of  cavity  radius  (Rc)  scaling  models  have  been  developed  since  the  beginning  of  the 
underground  nuclear  testing  era.  The  most  widely  used  by  the  seismology  community  are  the  ones  developed  by 
Denny  and  Johnson  (1991)  (DJ91)  and  by  Heard  and  Ackerman  (1967)  (HA67).  At  nominal  depth  of  burials  (DoBs) 
the  predictions  of  these  two  scaling  relations  are  in  quite  fair  agreement  between  them  and  also  with  experimental 
data.  However  when  the  DoB  is  increased  beyond  the  nominal  values,  the  predictions  steadily  diverge  and  differ  by 
51%  at  1000  scaled  meters  for  a  granite  medium.  As  a  result,  yield  estimates  made  for  the  25  May  2009  North 
Korean  test  based  on  seismic  magnitudes  depend  strongly  on  the  chosen  scaling  model.  Hydrodynamic  calculations 
involving  realistic  material  response  models  offer  a  viable  alternative  to  characterize  Rc  scaling  for  a  range  of  scaled 
DoB  (SDoB)  where  limited  data  from  past  nuclear  tests  exist.  Results  of  such  calculations  are  presented  for  a  granite 
medium  with  a  material  response  validated  by  modeling  four  phenomenological  criteria  for  past  nuclear  tests  in 
granite  (free  field  velocity  waveforms,  energy  partitioning  into  the  seismic  wavefield,  velocity  attenuation  profiles, 
and  measured  Rc).  The  numerical  results  presented  in  this  paper  unambiguously  favor  the  Rc  scaling  model  of  DJ91. 
Lower  bounds  on  yield  and  DoB  of  the  North  Korean  test  are  constrained  by  predictions  of  an  SDoB  threshold  for 
free  surface  damage  from  2-D  simulations  since  no  such  reported  damage  was  observed  for  this  test.  The  D J9 1 
model  indicates  the  minimum  yield  and  DoB  for  the  25  May  2009  North  Korean  test  is  5.7  kilotons  and  375  m. 

Background.  Traditional  approaches  to  seismic  yield  estimation  depend  heavily  on  empiricism  and  on  calibration  of 
path  effects  for  test  sites  around  the  world.  Even  the  best  efforts  to  calibrate  test  sites  were  not  enough  to  monitor 
compliance  with  the  Threshold  Test  Ban  Treaty  (TTBT),  and  treaty  protocols  had  to  be  re-negotiated  to  allow  for 
Joint  Verification  Experiments  (JVEs)  on  both  U.S.  and  Soviet  test  sites  (Douglas  and  Marshall,  1996).  For  broad 
area  monitoring  purposes,  yield  estimation  is  further  challenged  because  empirical  approaches  are  plagued  by 
inadequate  calibration  data  and  by  uncertainties  about  the  source  emplacement  which  may  not  conform  to  standard 
nuclear  containment  practices  on  which  seismic-derived  empirical  methods  are  based. 

Results  of  the  JVEs  validated  seismologists'  claims  of  the  existence  of  a  test  site  body  wave  magnitude  (nib)  bias 
between  U.  S.  and  the  former  Soviet  Union  test  sites  in  Nevada  and  Semipalatinsk.  The  use  of  remotely-collected 
seismic  data  to  calibrate  test  sites  was  vindicated,  and  it  was  shown  that  calibration  issues  can  be  overcome  even  for 
broad-area  monitoring  applications.  However,  the  concerns  about  non-standard  containment  practices  persist  and 
need  to  be  addressed  as  experience  with  the  DPRK  tests  shows. 

A  plot  of  yield  vs.  DoB  tradeoff  curves  (TOCs)  for  the  25  May  2009  DPRK  test  is  shown  in  Figure  1.  The  TOCs  are 
based  on  empirical  Rc  scaling  relationships  developed  by  HA67,  used  for  the  Mueller  and  Murphy  (1971)  model, 
and  DJ91.  These  TOCs  were  obtained  from  an  analysis  of  mt  and  moment  magnitude  (Mv),  while  accounting  for 
possible  nib  bias  and  the  effects  of  source  medium  elastic  properties  at  the  DPRK  test  site  (Patton,  2011).  Both 
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scaling  relations  are  consistent  with  the  measurements.  The  ambiguity  over  Rc  scaling  and  the  burial  depth  for  this 
test  translates  into  growing  yield  uncertainty  as  DoB  increases. 


DoB  (m) 

Figure  1.  Yield  :  DoB  tradeoff  curves  for  the  North  Korean  2009  test  (after  Patton,  2011).  Tradeoff  curves 
based  on  HA67  and  DJ91  Rc  scaling  equations  are  shown  with  black  and  red  lines,  respectively. 


Non-linear  material  responses  and  strong-motion  phenomenology  of  the  source  medium  control  Rc  evolution  and 
play  a  significant  role  in  seismic-wave  generation  of  buried  explosions.  Nuclear  testing  was  dominated  by  standard 
containment  practices.  As  such,  much  better  agreement  is  seen  between  Rc  scaling  models  for  SDoB  consistent  with 
these  practices.  Due  to  limited  Rc  data  for  non-standard  burial  depths,  hydrodynamic  computations  involving 
realistic  material  response  models  offer  a  viable  alternative  to  answering  questions  about  Rc  evolution  and  statics.  In 
this  paper,  such  questions  are  investigated  through  a  series  of  computations  for  a  granite  medium.  The  results  will 
test  the  veracity  of  Rc  scaling  models  and  will  even  help  constrain  the  SDoB  of  the  2009  DPRK  test. 

Hydrodynamic  Modeling  in  Granite.  A  wide  range  of  physical  and  thermodynamic  phenomena  control  shock 
wave  propagation  for  underground  nuclear  explosions.  Four  different  phenomenological  zones  can  be  identified 
according  to  the  material  behavior.  The  first  zone,  usually  called  the  source  region,  is  where  the  rock  and  device 
material  are  vaporized  due  to  the  huge  amount  of  energy  released  from  the  explosion.  As  the  shock  wave  propagates 
outward,  the  energy  dissipates  due  to  mechanical  and  geometric  factors.  At  a  certain  range  the  amount  of  energy  is 
no  longer  sufficient  to  vaporize  the  rock  but  is  sufficient  to  melt  it  in  a  second  zone  called  the  cavity  region.  As  the 
shock  wave  propagates  further,  the  rock  material  no  longer  melts,  but  it  undergoes  significant  irreversible 
deformations  and  fracturing.  This  defines  a  third  zone  called  the  permanently  deformed  region.  The  material  starts  to 
behave  elastically  at  a  radius  roughly  an  order  of  magnitude  larger  that  R^  therefore  defining  the  fourth  zone  called 
the  elastic  region.  In  the  numerical  results  Rc  is  considered  to  be  the  radius  that  defines  the  boundary  between  the 
second  and  the  third  phenomenological  zones. 

To  properly  account  for  wave  propagation  across  all  the  zones  and  also  for  the  correct  energy  deposition  in  the 
elastic  zone,  a  series  of  key  design  constraints  are  identified.  These  constraints  are  obtained  from  the  available 
experimental  data  for  Hardhat  and  Piledriver  explosions  detonated  in  granite  on  Climax  Stock  of  the  Nevada  Test 
Site,  (Rimer  et  ah,  1990;  Murphy,  1978;  Antoun  et  ah,  2001;  Perret  and  Bass,  1974),  and  are  to  be  employed  as 
design  parameters  for  the  purpose  of  developing  an  appropriate  material  model.  The  constraints  selected  are:  (1) 
velocity  and  displacement  waveforms,  (2)  energy  partition  for  each  zone  defined  above,  (3)  peak  velocity 
attenuation  profiles  and  (4)  Rc. 

The  three  most  common  approaches  for  modeling  the  source  region  are:  iron  gas  model,  bubble  model  with  ideal 
gas,  and  bubble  model  with  tabulated  equation  of  state  (EOS).  The  bubble  model  considers  that  the  energy  of  the 
explosion  is  uniformly  distributed  over  a  vaporization  sphere  with  radius  in  meters,  Ry  -  2.0  W  Schroeder  (1974), 
where  W  is  the  yield  of  the  explosion  measured  in  kilotons.  The  iron  gas  model  assumes  that  after  the  explosion  an 
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iron  gas  is  formed  occupying  a  spherical  region  around  the  shot  point.  The  initial  energy  of  the  explosion  is 
uniformly  distributed  over  this  spherical  region,  which  has  a  radius  smaller  than  Ry.  Different  material  models  have 
been  used  in  the  past  to  describe  the  behavior  of  the  “bubble”  or  source  region.  The  two  most  common  material 
models  employed  for  this  purpose  are:  ideal  gas  with  a  density  equal  to  the  density  of  the  rock  under  consideration 
(Antoun  et  al.,  2001),  and  a  tabular  EOS,  i.e..  Sesame  tabular  EOS  (Abdallah  et  ah,  1980). 

The  three  approaches  for  modeling  the  source  region  were  tested  against  Hardhat  and  Piledriver  observations,  and 
none  satisfied  all  of  the  design  constraints.  The  results  for  the  bubble  model  with  an  ideal  gas  matched  constraints  #1 
and  #3,  but  over  predicted  #4  and  the  energy  distribution  #2  was  not  realistic.  After  the  shot,  a  very  small  amount  of 
the  total  energy  (around  4%)  was  in  the  source  region,  while  in  reality  at  least  30%  of  the  energy  released  is  needed 
just  to  vaporize  the  granite  contained  within  Ry.  The  iron  pill  and  bubble  models  with  the  Sesame  EOS  gave 
reasonably  good  results  for  #4  and  energy  distribution  #2,  but  both  under  predicted  #1  and  #3.  These  EOS  models  do 
not  fiilly  account  for  the  first  principle  physics  of  shock  wave  continuum  from  the  cavity  expansion  to  final  energy 
deposition. 

The  importance  of  correct  energy  deposition  and  cavity  dynamics  required  finding  a  different  EOS  model  for  granite 
media.  Tillotson  (1962)  introduced  a  dependency  on  the  specific  internal  energy  of  the  material  and  accounted  for 
phase  changes  from  solid  to  vapor.  A  series  of  parameters,  based  on  a  modified  version  due  to  Melosh  (1989),  are 
used  to  describe  the  behavior  of  the  material  across  the  different  phases.  The  source  modeling  was  then 
accomplished  by  coupling  the  Tillotson  EOS  with  the  bubble  model. 

The  material  model  is  completed  by  combining  the  EOS  with  an  appropriate  strength  model.  The  main  component 
of  the  strength  model  for  granite  is  the  “yield  surface”,  which  describes  how  the  shear  strength  changes  with 
confinement  pressure.  The  adopted  undamaged  yield  surface  (Yundamage)  was  derived  from  Fossum  and  Brannon 
(2004),  where  modifications  were  introduced  to  match  the  field  data  for  Hardhat  and  Piledriver. 

The  material  model  also  incorporates  a  shear  strain  softening  behavior,  which  follows  the  one  presented  by  Rubin  et 
al.  (2000).  The  undamaged  yield  surface  curve  is  scaled  down  as  a  function  of  the  strain  softening  factor  (fso/t)  as 
follows 

^  “  ^undamagef  soft  ’  f  soft  —  1  “  ^ 

V  y 

where  fsoft  is  the  softening  factor,  a  and  b  are  material  fitting  constants,  and  Q  =  min(^^  I shear  damage, 

where  Sp  is  the  plastic  shear  strain  and  is  the  value  of  the  plastic  shear  strain  for  which  the  damage  is  maximum, 
i.e.  Q  is  one  and  therefore  =  \  —  a  . 

Experimental  data  from  Hardhat  and  Piledriver  explosions  were  used  as  a  reference  for  adjusting  the  material  model 
parameters.  Much  work  has  been  done  in  the  past  to  simulate  these  two  explosions  with  different  material  models 
(Antoun  et  al.,  2001;  Stevens  et  al.,  2000;  Stevens  et  al.,  2001).  However,  to  the  authors’  knowledge,  constraint  #2  is 
yet  to  be  reported  in  this  kind  of  study.  In  the  following  section,  the  objective  of  satisfying  all  design  constraints  is 
accomplished  using  the  modeling  approach  just  described  for  1-D  simulations. 

Fits  to  freefield  ground  motion,  energy  distribution  and  Rc.  Figure  2  shows  a  comparison  of  the  calculated 
against  the  observed  free  field  velocity  waveforms  for  the  Hardhat  event  as  a  function  of  time. 
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Hardhat-  Free  Field  Velocity  Profile 

Range  =  184  m 


Hardhat  -  Free  Field  Velocity  Profile 

Range  =  239  nn 


Time  [s] 


Time  [s] 


Figure  2.  Material  model  validation:  HARDHAT:  Comparison  of  calculated  vs.  observed  velocity  profiles  for 
different  ranges.  The  calculated  velocity  profiles  were  obtained  using  the  material  model  with  just 
strain  softening  behavior. 


At  the  end  of  the  Hardhat  simulation  the  energy  distribution  per  zone  was:  zone  1:  46%,  zones  2  and  3:  50%,  and 
zone  4:  4%.  The  amount  of  energy  deposited  in  zone  4  gives  the  seismic  efficiency  of  the  underground  nuclear 
explosion  which  compares  favorably  with  1.8%  for  the  granite  shot  Shoal  reported  by  Mueller  (1969).  This  high 
energy  value  obtained  from  the  Hardhat  simulation  and  the  slight  over-estimation  of  peak  velocities  on  Figure  2  can 
be  explained  by  the  fact  that  many  dissipative  factors  are  not  currently  considered  in  the  computational  model,  such 
as  pre-existing  faults  and  fractures. 

As  the  energy  propagates  across  the  computational  model  a  certain  amount  of  work  due  to  compression  is  done  on 
the  cells  of  the  model.  This  work,  which  is  essentially  given  by  p  dv,  where  dv  is  the  volume  change,  increases  the 
specific  internal  energy  of  the  cell,  therefore  increasing  its  temperature.  When  the  specific  internal  energy  reaches  a 
given  value,  the  material  starts  to  melt.  The  position  of  the  last  melted  cell  from  the  center  of  the  source  was  the 
criterion  used  to  identify  where  the  boundary  of  the  cavity  area  is  located.  The  calculated  and  measured  values  of  Rc 
for  Piledriver  are  39.4  and  44.5  respectively  and  for  Hardhat  are  17.9  and  19.4  m.  Differences  between  calculated 
and  measured  Rc  of  only  1 1.5%  and  7.7%  respectively  are  indication  of  fairly  good  agreement.  It  is  worth  noting  that 
the  predictions  of  HA67  and  DJ91  for  Hardhat  are  21.4  and  18.4  m  respectively. 

Effects  of  Scaled  Depth  of  Burial.  In  this  section,  the  effects  of  over-burden  on  Rc  evolution  and  on  material 
damage  for  a  2-D  axially  symmetric  computational  model  are  investigated  utilizing  the  material  models  discussed 
above.  A  1-kt  source  is  placed  at  different  DoB  (125,  250,  500  and  1000  m)  in  a  homogeneous  granite  medium  with 
Young  modulus:  70.22  GPa,  density:  2680  kg/m^  and  Poisson  ratio:  0.226  (P-wave  velocity:  5495  m/s,  S-wave 
velocity:  3269  m/s).  The  free  surface  is  flat,  i.e.  no  topography.  The  influence  of  the  gravitational  field  is  taken  into 
account  in  the  model,  i.e.  a  lithostatic  stress  state  is  prescribed  as  an  initial  condition. 

The  final  value  of  Rc  is  recorded  for  each  DoB  and  compared  to  the  predictions  of  scaling  relationships  due  to  HA67 
and  DJ91  as  shown  in  Figure  3.  Plotted  is  HA67’s  empirical  scaling  equation  for  R^  the  same  as  equation  (19)  in 
Mueller  and  Murphy  (1971),  and  equation  (39)  of  DJ91.  The  scaling  equations  show  good  agreement  for  SDoB  near 
120  m/kt^^^.  Most  Rc  measurements  on  which  the  empirical  scaling  equations  were  based  are  for  explosions 
detonated  at  these  nominal  containment  depths.  Far  fewer  measurements  were  available  for  large  SDoB  where  the 
predictions  steadily  diverge  to  the  point  that  they  differ  by  51%  (1 1.2  versus  7.4  m)  at  1000  m/kt^^^.  The  numerical 
results  are  in  much  better  agreement  with  DJ9rs  scaling  equation,  even  for  quite  deeply  buried  events. 

Another  effect  of  SDoB  is  the  amount  of  permanent  damage  due  to  shock  wave  interactions  with  the  free  surface. 
Figure  4  shows  shear  damage  profiles  for  different  SDoBs.  At  a  nominal  SDoB,  there  is  a  substantial  amount  of 
shear  damage  on  the  free  surface,  while  damage  decreases  as  SDoB  increases,  as  expected.  It  can  be  inferred  that  a 
SDoB  of  209  m/kt^^^  is  the  limit  where  almost  no  shear  damage  is  noticeable  for  a  uniform  granite  medium  with 
material  properties  specified  above  and  a  flat  free  surface.  It  is  worth  noting  that  extensive  free  surface  damage  (up 
to  a  range  of  753  m)  in  the  form  of  spallation  was  observed  for  the  Piledriver  event,  (Bache,  1976). 
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Figure  3.  Cavity  radius,  Rc,  as  a  function  of  SDoB  for  a  1  kt  explosion  in  granite. 


Figure  4.  Shear  damage  for  different  SDoBs.  Yield  =  1.0  kt. 


Material  Model  Enhancements.  As  part  of  the  ongoing  research  effort  the  material  model  developed  last  year  was 
enhanced  by  incorporating  shear  hardening  on  top  of  the  already  implemented  shear  softening  behavior.  The 
hardening  process  adopted  in  this  work  is  described  in  detail  by  Rubin  et  al.  (2000).  The  working  yield  surface  is 
then  obtained  as  follows 


7  = 


V  f  f 

undamage  J  soft  J  hard 


(2) 


where  Y undamage  IS  thc  yield  surface  of  the  intact  material,  is  given  by  equation  (1)  dtndfhard  is  the  hardening  factor 
which  is  given  by 


/w  =c  +  (l-c)exp(-f/) 


(3) 


where  c  is  the  amount  of  hardening  for  the  material,  Sp  is  the  equivalent  plastic  strain  and  J  is  a  parameter  that 
defines  how  fast  the  shear  hardening  is  applied  as  a  function  of  Sp.  With  the  introduction  of  the  hardening  factor  the 
response  of  the  material  as  it  undergoes  plastic  deformation  in  shear  can  be  divided  into  two  phases:  in  the  first 
phase  the  material  hardens  until  it  reaches  maximum  shear  strength.  As  the  equivalent  plastic  strain  grows  softening 
behavior  starts  to  be  more  preponderant  than  the  hardening  and  therefore  the  shear  strength  degrades  up  to  a  defined 
minimum  as  shown  in  Figure  5a.  In  Figure  5b  the  corresponding  yield  surfaces  are  shown. 
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a) 


b) 


Figure  5.  a)  Softening  and  hardening  factors  as  a  function  of  the  equivalent  plastic  strain,  b)  Yield  surfaces 
utilized  in  the  improved  material  model. 


In  order  to  adjust  this  newly  developed  material  model  a  set  of  numerical  simulations  were  conducted  for  the 
Hardhat  event.  The  results  obtained  for  the  radial  velocity  waveforms  at  ranges  184  m  and  239  m  are  shown  in 
Figure  6.  The  matching  between  the  calculated  results  and  the  observed  values  is  excellent  and  the  improvement  can 
be  readily  appreciated  (see  Figure  2).  The  final  cavity  radius  obtained  from  the  simulations  was  19.0  m,  which  is  in 
much  closer  agreement  (2%  difference  with  the  measured  value,  19.4  m)  than  the  results  obtained  using  the  material 
model  developed  last  year,  17.9  m. 


Hardhat  -  Free  Field  Velocity  Profile 

Range  =  184  m 


Time  [s] 


Hardhat  -  Free  Field  Velocity  Profile 

Range  =  239  m 


Time  [s] 


Figure  6.  HARDHAT:  Comparison  of  calculated  vs.  observed  velocity  profiles  for  different  ranges.  The 

calculated  velocity  profiles  were  obtained  using  the  improved  material  model  (with  strain  hardening 
and  softening  behavior). 


CONCLUSIONS  AND  RECOMMENDATIONS 


A  lower  yield  limit  for  the  2009  North  Korean  test  can  be  established  by  combining  the  above  results  with  the  TOCs 
discussed  in  earlier  in  the  paper.  Satellite  images  of  the  2006  and  2009  tests  showed  no  evidence  of  free  surface 
disturbances  caused  by  spallation  (Schlittenhardt  et  al.,  2010;  Begnaud  et  al.,  201 1).  Figure  7  shows  curves  of 
constant  SDoB  corresponding  to  125,  209  and  250  m/kt^^^  plotted  with  the  TOCs.  The  curve  for  209  m/kt^^^  is  the 
calculated  estimate  of  the  threshold  for  shear  damage  expected  at  the  free  surface  for  an  explosion  in  granite:  i.e.,  if 
the  explosion  was  buried  much  shallower  than  209  m/kt^^^,  the  effects  of  spallation  and  free  surface  damage  should 
have  been  great  enough  to  be  noticed  from  satellite  imagery.  This  constraint  means  that  a  combination  of  DoB  and 
yield  should  be  located  to  the  right  of  this  curve.  Since  the  numerical  results  favor  the  Rc  scaling  equation  of  DJ91 
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over  HA67  (Figure  3),  estimates  of  the  minimum  DoB  and  yield  for  the  2009  DPRK  test  are  375  m  and  5.7  kt, 
respectively.  Greater  DoB  and  larger  yields  should  follow  the  trajectory  defined  by  the  DJ91  TOC  curve. 


DoB  (m) 


Figure  7.  Yield  vs.  DoB  tradeoff  curves  for  the  North  Korean  2009  test  (after  Patton,  2011).  Tradeoff  curves 
based  on  HA67  and  DJ91  Rc  scaling  equations  are  shown  with  black  and  red  lines,  respectively. 
Trajectories  of  constant  SDoB  are  shown  by  yellow,  green,  and  brown  lines. 


In  summary,  hydrodynamic  simulations  play  an  important  role  in  investigating  the  effects  of  non-linear  material 
response  and  strong-motion  phenomenology  on  Rc  scaling  and  seismic-wave  generation  for  source  emplacements 
with  limited  experience  or  data  from  past  nuclear  tests.  The  focus  of  this  study  was  to  predict  final  Rc  for  explosions 
detonated  in  a  granite  medium  and  to  determine  a  threshold  for  permanent  damage  on  the  free  surface  as  SDoB 
increases.  The  results  agree  with  the  empirical  Rc  scaling  model  of  DJ91  and,  along  with  a  previously  developed 
yield  :  DoB  tradeoff  curve,  constrain  lower  bounds  on  the  yield  and  DoB  of  the  25  May  2009  DPRK  test. 

3D  Topographic  Effects.  Refined  3D  simulations  including  surface  topography,  complex  structures  and 
morphology  will  be  conducted  to  augment  the  understanding  of  Rc  scaling  and  seismic-wave  generation  for  realistic 
explosion  source  regions.  The  first  stage  of  this  next  step  will  be  to  develop  a  computational  model  taking  into 
account  the  topography  of  the  Nevada  National  Security  Site  (NNSS)  area  15-b  (Figure  8)  and  to  investigate  the 
effects  of  introducing  topography  on  the  generation  of  shear  waves.  Preliminary  results  will  be  presented  at  the 
conference. 


Figure  8.  Left:  Digitized  Elevation  Map  of  the  NNSS  Area  U-15b.  Right:  Finite  volume  model  with 
topography. 
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Near  to  Far  Field  Propagation.  As  part  of  this  year’s  efforts  we  have  begun  the  process  of  coupling  the 
hydrodynamics  output  to  a  far  field  seismic  wavefield  propagation  code  SPECFEM3D.  SPECFEM3D  is  a  spectral 
element  code  that  allows  for  complex  topography,  full  othorhombic  anisotropy  and  intrinsic  attenuation 
(Komatitsch,  1997;  Komatitsch  and  Vilotte,  1998).  Our  current  model  has  the  topography  of  Area  15  at  the  Nevada 
National  Security  Site  and  a  two  layer  model  of  the  weathered  and  intact  Climax  Stock  granodiorite.  Careful 
analysis  of  the  shear  wave  field  generated  near  source  will  be  compared  with  that  generated  in  the  far  field. 

The  above  analysis  can  then  be  followed  with  the  development  of  a  computational  model  that  considers  the  real 
topography  of  the  North  Korean  test  site  (Figure  9).  In  so  doing,  further  refinement  of  the  minimum  yield  and  depth 
of  burial  predicted  earlier  can  be  achieved. 
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Figure  9.  A  photograph  of  the  location  of  the  East  and  West  portals  at  the  North  Korean  test  site,  Pung’gye,  is 
shown. 
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